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Abstract 

The main purpose of this work is to simulate two-phase flow in the form of immiscible displacement through 
anisotropic, three-dimensional (3D) discrete fracture networks (DFN). The considered DFNs are artificially gener- 
ated, based on a general distribution function or are conditioned on measured data from deep geological investigations. 
We introduce several modifications to the invasion percolation (MIP) to incorporate fracture inclinations, intersection 
lines, as well as the hydraulic path length inside the fractures. Additionally a trapping algorithm is implemented 
that forbids any advance of the invading fluid into a region, where the defending fluid is completely encircled by the 
invader and has no escape route. We study invasion, saturation, and flow through artificial fracture networks, with 
varying anisotropy and size and finally compare our findings to well studied, conditioned fracture networks. 
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1. Introduction 

Fracture networks are abundant in natural geological formations. Characterizing and investigating them is among 
the most challenging problems in geosciences, ranging from ground-water hydrology, subsurface oil, isolation of 
nuclear and hazardous wastes or captured CO2, to geotechnical problems involved in tunneling and cavern excavations. 
The identification and assessment of fractures as capillary barriers, the understanding of flow and transport in fracture 
networks, as well as the development of mathematical models is crucial for predicting the hydraulic behavior of 
fractured geological systems. Due to the geometrical complexity of the natural system, analytical solutions are not 
available, but also it is not feasible to describe the systems in detail |nil2|. To be able to study natural systems, their 
structure, as well as the occurring processes, they have to be represented by conceptual models. Numerical solutions 
have to be applied which involve spatial and temporal discretization of the problem and usage of efficient and stable 
algorithms. 

One approach to modeling flow in fractured rock is to represent the fracture network within the rock as a three- 
dimensional discrete fracture network (DFN). The DFNs are either given on the basis of core samples and measure- 
ments on exposures (tunnel walls or outcrops), or they emerge from planar, irregular, ^-polygonal fractures that are 
oriented around the mean direction, which may follow for example a Fisher distribution O. The latter are considered 
to be a solid representation of the ones experimentally measured |4-7 |. One possible approach for modeling the dis- 
placement of immiscible fluids like gas, oil, or water inside DFN models, is the so-called invasion percolation (IP). 
IP is a dynamic growth process that generates phase structures through a set of rules, which embody the physics of 
immiscible displacement of fluids within a random field, such as a network of pores or fractures. IP was proposed 
by Lenormand and Bories \S \ and Chandler et al |9 | in the early 1980s, as well as Wilkinson and Willemsen who 
coined the name IP ifTOl . They emphasized that IP is a modified form of ordinary percolation (OP) ifTTIl . only with a 
well-defined sequence of invasion events. In IP it is assumed that capillary forces dominate the viscous ones. In other 
words, both gravity and viscous eff'ects are neglected. Nevertheless, IP is a valid approximation for the flow physics 
and well suited to describe the slow immiscible displacement of two fluids in porous media or fracture networks. 
In this work, physical modifications to the model are proposed to study the gas-water displacement process in the 
DFNs. IP has been modified (MIP) in many ways in the past to incorporate more geometrical and physical details of 
the network like force fields such as gravity |[T2lfT5ll , centrifugal forces O, viscous forces ifTTl . as well as capillary 
smoothing mechanisms in porous media (131 llSl and inside single fractures {191 . In this work, we have modified the 
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IP model to a more realistic description of the invading progress. To decide if the invader will penetrate the fracture 
adjacent to the invading front, not only the action of a fracture as a capillary barrier, but also the contact profile to 
this fracture and its inclination with respect to the flow direction, as well as the hydraulic path length from entry to 
exit will be considered. As the invasion proceeds, part of the invasion percolation cluster boundary might completely 
encircle a region of the fracture network. This phenomenon, named trapping, calls for a sophisticated IP that forbids 
any advance into this area. The IP or MIP routine is appropriately called trapping IP 1 10] or trapping MIP. To study 
leakage rates and pressure fields above percolation, after pruning dead ends (backbone of the percolating cluster), the 
cubic law for the remaining fracture network is solved in the form of a linear hydraulic system. The global system 
matrix is composed of invaded fractures as elements and the connections as nodes using a simple channel model 
following Cacas et al. 1201 . 

The paper is organized as follows: Sec. [2] provides information on 3D fracture networks in general and describes 
the generation of discrete fracture network (DFN) models based on either a general distribution function or on previ- 
ously measured distributions. The physical modifications of the invasion percolation (MIP) procedure including the 
trapping phenomenon are introduced in Sec. [3] followed by the explanation of the typical invasion scenarios in Sec.|4] 
Sec.|5]is devoted to the physics of flow in single fractures and subsequently in fracture networks. In Sec.[6]the results 
of simulations in artificially generated fracture networks with Fisher distributed fracture orientations are presented and 
complemented by a size dependence study of the relevant properties in Sec. [7] Finally, results are compared with those 
from a model conditioned on data from geological investigations |21| in Sec. [8] The paper closes with conclusions 
and an outlook in Sec.[9l 

2. Fracture networks 

Natural geological systems are highly complex in terms of the geological structure as well as the physical processes 
occurring within them, making analytical solutions unfeasible. Hence, a simplifying transformation of the natural 
system to a numerical model is required |2|. The natural fracture system has to be transformed to a discrete fracture 
network (DFN) that is composed of individual, intersecting fractures with similar hydraulic characteristics as the real 
system. To be able to make systematic studies, that are comparable to the state of the art knowledge, the artificial 
DFNs are constructed in a similar manner to those of Refs. |22, 23 1. 

The hydraulic properties of fracture network systems are mostly determined by the permeability of the fractures 
and matrix. In our conceptual model, the permeability is exclusively determined by fractures. The matrix between the 
fractures is not considered, since its permeability compared to the ones of the fractures is assumed to be extremely low 
(see | 24 |). The hydraulic properties of the DFN are typically characterized by a fracture size, fracture permeability, 
fracture orientation, and fracture density distribution. Individual fractures can be seen as a two-dimensional porous 
medium with flow inside the fracture plane. Variations of pressure or velocity across the height /z, called the aperture 
of the fracture, can be averaged out. The aperture is much smaller than the other two dimensions of the fracture and 
has an essential influence on the flow and transport processes in fracture networks. In the DFN model, the aperture is 
described by the parallel-plate concept (Fig. [T]), that is the simplest model of flow through a single fracture with the 
possibility for exact solutions of the hydraulic conductivity. 
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Figure 1: Parallel-plate concept: a) From natural single fracture to local parallel plates, b) From local parallel to one single plate with mean aperture 
h (distance between the plates) following SiiB |2 1. 

2.1. Artificial fracture network models 

The artificial DFN is composed of Np planar, regular or irregular, convex polygons with Ny vertices. Each fracture 
is approximated by an individual polygon / (with / = 1 . . .Np) inscribed in a disk of radius Ri. This approximation 
comes from the parallel plate model, assuming that rigid walls are smooth and plain. The radii Ri of the circumscribed 
circles are randomly chosen in the interval [^min,^max], and the number of vertices Ny is picked in the interval 
[3, A^max], both from a uniform distribution. The vertices are distributed on the perimeter of the circle with uniformly 
distributed angles Ofi, 0^2, . . . , aj^^, all in the interval [0, In]. A simplification can be made in terms of choosing regular, 
equal sized (monodispersed) polygons on circles with fixed radii Ri = 7?max, as well as an identical number of vertices 
Ny and subsequently equal angles aj = InlNy. The aperture is in this case just the ratio of the pore volume Vpoiy,/ over 
the area Apoiy,/. The pore volume Vpoiy,/ is calculated using Vpoiy,/ = A^^^^ Vchar, with the uniformly distributed random 
variable Zi e [0, 1] and the characteristic pore volume Vchar, which is a given constant depending on the fracture system 
under consideration. Detailed investigations of fractures suggest that aperture varies within each fracture (e.g. f25l). 
However, within the current model we assume the fracture aperture to be uniform over each polygon. Hydrogeological 
investigations, where transmissivity or equivalently hydraulic aperture is assumed to be uniform within fractures and 
is determined by pumping tests, typically show a long-tailed aperture distribution for all fractures, often represented 
as a log-normal or Pareto distribution |26|. 

Once all polygonal fractures are generated, their centers are placed in the simulation box of size L^, using a 
uniform distribution before orientations are adjusted. The homogenization scale L is chosen such that it significantly 
exceeds the size of the largest possible fracture ^ by L » Jy^, with dk = <imax the diameter of the circumscribed circle 
of fracture k. The fracture orientation is represented by the normal vector of the fracture plane. The distribution of 
a random set of normal vectors can be described by means of a probability density function f(0, cp), expressed in 
standard spherical coordinates 6 and 0. The simplest probability density function would be the uniform distribution 
of the normal vectors on the unit sphere. However, in order to study anisotropic, three-dimensional (3D) DFN, which 
are closer to reality, the Fisher distribution can be used ||23]|27l. This distribution is the analogue of the Gaussian 
distribution on the sphere 1 3 1 and can be deduced in polar coordinates e [0, tt] and g [0, 27i] as follows: 

f{6, 0, a:) = : sin 6 exp (k [cos ^0 cos 6 + sin ^0 cos(0 - 0o)l) , (1) 

An smh k 

with /c > denoting the concentration or dispersion parameter. For the particular case where the initial coordinates 
^0 and 00 are equal to zero, the Fisher distribution reduces to 

f{6, k) = - — ^- — sin 6 exp (a: cos 6) . (2) 
4n sinh k 

This distribution is rotationally symmetric around the initial mean direction, which coincides with the z-axis. The 
greater the value of k, the higher the concentration of the distribution around this axis. The distribution is unimodal 
for /c > and is uniform on the sphere for /c = 0. 
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Examples of generated artificial 3D DFNs with diff'erent k values are given in Figs.|2]a) and b). Note that only when 
fractures are highly connected, the system behaves like a continuous medium. Since the interest lies in the movement 
of a fluid phase in the network, the fracture density pfr of a DFN, being the number of fractures per volume, has to 
be high enough to build a spanning cluster over the entire domain. The connectivity of the DFN can be resembled 
by the dimensionless density or the so-called concentration p . Huseby et al. 1 22 1 have shown, that the concentration 
p equals the average number of intersections per fracture, p' = ^Mnt). The concentration is thereby a decreasing 
function of k 1231 . Hence, a much smaller average number of intersections per fracture is expected for the highly 
anisotropic DFN (see Fig.|2]b)) than for the nearly isotropic one (see Fig.[2]a)). Furthermore, Khamforoush et al. | 27 1 
noted that the asymptotic values (for the infinite system L oo) of percolation thresholds (of OP) for anisotropic 
DFNs {pc = p'c) are in 

2.1 < Pc,oo,.,A' Pc,oo,y,A ^ 2.3 and 2.3 < p^^^^,, < 2.44 . (3) 

p^^ denote critical concentrations in the range of < a: < 50 with the mean direction of the Fisher distribution chosen 
to be the z-axis. 

2.2. Realistic fracture network models 

In addition to the artificial networks described in the previous section it was thought useful to consider more 
realistic models conditioned on field data. Within the study, network models of the Excavation Damage Zone around 
tunnels in Opalinus Clay (not described here) and a large-scale model of the fracture system within a fractured marl 
were considered. The large-scale model is a simplified version of that described in Mazurek et al | 21 1 and Nagra 
[281 . It was selected as a model likely to show a well-connected 3D fracture system suitable for comparison with 
the analytic models. The model is conditioned on field data acquired from deep boreholes in the Valanginian Marl at 
Wellenberg in Central Switzerland during Nagra's investigations there in the early 1990s The model includes 
fracture "sets" corresponding to the inventory of water conducting features (WCFs) identified in the boreholes. The 
WCFs included cataclastic faults zones over a wide range of length scales and smaller structures including thin discrete 
shear zones, joints and isolated limestone beds. The full model also included a small number of larger-scale limestone 
banks that had been identified but these were not in the stochastic model used here. The length distribution of the 
cataclastic zones was based on the observed thickness and roughly followed a power-law form with exponent between 
2.3 and 2.6. Feature orientation was based on observations from core and showed a strong relationship to the folding 
axis of the rock but includes discordant sets resulting in a well connected network. Hydraulic aperture distributions 
were calculated from the observed transmissivity measured by pumping tests and fluid logging analyses. 

The Wellenberg DFN model shown in Fig.|2]c) is comparable to the 3D randomly generated artificial fracture 
networks and is built similarly in the form of a system with L = 100 m. Large natural faults, spanning over hundreds 
of meters, are thereby tessellated into smaller fractures, represented by polygons with side length of 7?/ = 10 m, using 
again the parallel plate model. The concentration in terms of the average number of intersections per fracture is 
p 15 and thus much higher than for the artificially generated DFN. This leads to a higher connectivity of the DFN 
and hence only a small amount of water is assumed to be trapped in the invasion percolation simulation afterwards. 
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Figure 2: a) Nearly isotropic (k = I) artificial 3D DFN with system size L = 10 dmax, where Jmax is the diameter of the circumscribed circle of 
the largest fracture in the DFN, containing 4000 fractures, b) Highly anisotropic (k = 50) artificial 3D DFN with the same size and number of 
fractures, c) Sample of the Wellenberg DFN with system size L = 100 m containing 12684 fractures. 



3. Physical modifications for IP 

IP in its simplest form fully invades an entire fracture, once the aperture threshold is reached and at least one 
intersecting fracture is filled. For fluid invasion, previous experiments have shown that the detailed behavior at fracture 
intersections are essential |29|. We use a modification to explicitly represent intersections and to implement rules at 
fracture intersections to mimic capillary barrier type behavior for the invading gas in form of an adjusted entry aperture 
h^^. . Given that some fractures of the fracture network models can be much bigger than others and that maybe only 
a small portion of them will be invaded, an additional modification considering the path length inside the fracture for 
the hydraulic aperture of a fracture is suggested. 

Consider an invaded fracture / from the invasion front, connected to a non-invaded fracture j, that is a feasible 
candidate for invasion. The resistance of fracture j is given by its entry aperture h^, which has to be at least as high 
as the current threshold value to allow for invasion. In the following, the aperture value of fracture j will be adjusted 
with respect to the eff'ects mentioned above: 

• The length of the line of intersection dc\ between fracture / and 7, which results from the connectivity calculation, 
is divided by the diameter of the circumscribed circle dj taken as the maximal width of fracture j (see Fig. [3]). 
Hence the discriminating factor concerning the line of intersection is 

nc = ^. (4) 

dj 

The corresponding factor for the invasion of entry fractures at the entry zone of the fluid is = 1 • 

• To consider the inclination of the fractures with respect to the flow direction, the aperture is corrected with a 
discriminating factor cos (0) (where < < 7r/2) as proposed by Sarkar et al. 1291 . The angle cp is either the 
contact angle between fracture / and j at the invading front (see Fig.jSj or the angle between the global inflow 
direction and the inclination of the considered entry fracture at the entry zone. This results for both cases in the 
discriminating factor of 

ni = cos (Ci0) , (5) 

with the adjustment parameter Q between and 1, to adjust the strength of the inclination eff'ect. Since cos (0) = 
for = 7r/2, the discriminating factor vanishes for fractures perpendicular to each other and flow between 
them would be prohibited. This is anticipated by setting Ci < 1. In the other extreme, no inclination adjustment 
is obtained by setting Cj = 0. 
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• The path length adjustment is done using the path length discriminating factor 

di 



np 



(6) 



where Jpi is the path length in fracture / from prior (invaded) entry node to the exit node to fracture j. di is the 
diameter of the circumscribed circle of fracture / being the maximal possible path length (see Fig.|3]). 



Finally the adjusted entry and adjusted hydraulic aperture of fracture j for invasion from fracture / is 

^adj,j = ncni/if and h\.=n^hf. 



(J) 



The ordinary IP is retrieved by the special case where np = nc = nj = 1 . The resulting total aperture is calculated 
using the rule for computing the equivalent aperture for fractures connected in series 1 30] 



1 



(8) 




Figure 3: Modified invasion percolation (MIP) with invaded fracture i at the invading front connected to the non-invaded fracture j (candidate for 
invasion). is the path length in fracture i from prior (invaded) entry node to the exit node to fracture j and di is the diameter of the circumscribed 
circle of fracture i. dd is the length of the line of intersection between fracture i and j, dj is the diameter of the circumscribed circle of fracture j 
and is the contact angle between the two fractures. 



3. 1 . Searching for traps 

As the MIP proceeds, part of the cluster boundary might completely encircle a region of the fracture network 
containing the defending fluid (see Fig. [4]). Given that the fluid is nearly incompressible, the entrapped region cannot 
be invaded. Hence the MIP version forbidding any advance into such an area is appropriately called trapping MIP 
(TMIP). Searching for traps and removing their perimeters from the list of potential invasion sites is cumbersome, but 
as a matter of fact one of the most important parts of IP simulations fST]. In early attempts, the entire lattice has been 
scanned using a Hoshen-Koplemann (HK) algorithm |31 1, which finds and labels all connected, water-filled regions 
that are disconnected from the outlet and scales as 0{N). Nevertheless, the computation time for the whole trapping 
algorithm is costly in computing time and scales as 0(N^) with the time for each lattice realization. For this reason 
Sheppard et al. | 32| designed a highly efficient algorithm whose execution time scales as 0(N In N). After insuring 
that trapping can occur, they used several simultaneous so-called breadth-first searches to update the cluster labeling. 
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This is similar to the strategy considered here, where we use the burning algorithm, first proposed by Herrmann et 
al. (331. To check the possibility of trapping, the last invaded fractures are included in a list of potentially trapped 
fractures. For each fracture in this list, a so-called depth-first search is used to check, for connectivity to the outlet. 
Fractures that are designated as trapped are labeled and removed from the invasion list. The efficiency of this algorithm 
stems from the fact that the number of traps decreases as a power law with the trap size | 31 1 and small traps can be 
detected quickly through local searches. 



4. Typical invasion scenarios 



Before the actual MIP simulation can be started, one has to determine, if the network percolates. Therefore a 
depth-first search called forward percolation is performed to ensure that the inlet is connected to the outlet. The inlet 
is defined by an entry zone and the outlet by an exit zone with a predefined volume expansion in 3D or width in 2D. 
Fractures located at least partially in one of these zones are labeled as entry or exit fractures, respectively. 

After the forward percolation, the trapping MIP is applied, till the breakthrough step is reached. At breakthrough, 
the largest interconnected cluster connects in- and outlet. This breakthrough state is summarized schematically in 
Fig. [4] with the finite clusters drawn by dashed lines and the infinite or percolation cluster by solid ones, connecting 
the entry and exit zone marked gray. The backbone is drawn in solid lines with links (chains of connected fractures) 
as lines, nodes (crossing points of the links) as hexagons and blobs (dense regions with more than one connection 
between two points like cycles or loops) as circles. The hatched area is representing the region containing trapped 
fractures. 




Figure 4: Two dimensional scheme of the system at the breakthrough representing the infinite and some finite clusters. The infinite cluster consists 
of a backbone (solid lines) and several dead ends (dotted line). Finite clusters are represented by dashed lines and trapped regions are hatched. 



Parts of the infinite cluster with dotted lines, connected to the cluster by only one bond, resemble dead ends. Since 
they do not contribute to the actual fluid (e.g. gas) flow after breakthrough, they can be pruned. Usually, the majority 
of the fractures form dead ends, but to simplify the scheme, only very few of them are shown in Fig. |4] The finite 
clusters, which are not part of the infinite cluster and will also not contribute to the flow process, are shown in dashed 
lines and can be pruned as well. What remains is the so-called backbone, being the portion of the infinite cluster 
that actually contributes to the fluid flow. In studies of gas leakage and underground waste repositories, the backbone 
is also called flow fracture network or leakage path. Some of the links are designated red or hot bonds, since no 
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alternate path exists and the spanning cluster would become finite, if they were cut. All definitions are summarized in 
the "links-nodes-blobs" picture similar to Fig.|4j introduced by Stanley |34|. 

The backbone is identified using again a burning algorithm but starting from invaded fractures at the exit zone and 
processing till it reaches invaded fractures in the entry zone. Since this process seems to run backwards in terms of 
the spatial direction given by the gas-water displacement simulated with MIP, it is called backward percolation. After 
this procedure, all fractures, that are not belonging to the identified infinite cluster are pruned. Also fractures with 
only one connection, which are neither inlet nor outlet fractures, are deleted. This is done recursively, till only the 
backbone remains. 

After reaching the breakthrough state, there are several scenarios, which can be considered in addition. Instead 
of stopping the MIP simulation at breakthrough, being followed by the backward percolation as discussed above, the 
invading process can continue, i.e. by assuming an infinite gas production rate at the inlet. The simulation can then 
be stopped either after /w// invasion of the system, meaning that no more water can be pressed out of the system and a 
specific water saturation due to trapping has been reached, or after reaching a certain threshold value, e.g. representing 
the highest possible input pressure at the inlet. Additional extensions could be pressure drops after breakthrough or 
above the percolation threshold after reaching again a certain threshold, to pretend a finite and maybe also dynamically 
varying gas production rate. In this work, the full invasion scenario is investigated (Secs.[6j [8]), additional to the usual 
simulation studies till breakthrough. All MIP variants or scenarios are followed by calculating post-breakthrough or 
collapsed-network properties such as the invaded volume, gas-water saturation relationship, pressure distribution, and 
flow rates. 



5. Estimates for permeability after breakthrough 

An idealized fracture network can be regarded as a network of tiny connected tubes. For small Reynolds numbers 
(of the order of one), Darcy's law can be applied, which relates mass flow to pressure gradient at the continuum scale. 
Darcy's law states, that the volumetric flow rate q of an incompressible fluid through a specimen of fractured porous 
material is given by 

.= ^, (9) 

/ll 

where A/? is the hydrostatic pressure difl'erence across the specimen, / is the length of the specimen, ji the dynamic 
viscosity, and A the cross sectional area. The constant of proportionality k is called the Darcy permeability of the 
material. Based on this simple relationship, mass balance leads to an equation for pressure in the domain of interest. 
The flow calculation through a single fracture using the parallel-plate concept from above yields the cubic law, namely 

wh^Ap 

where q again denotes the volumetric flow rate, Ap the pressure difl'erence, jd the dynamic viscosity, w the width of 
the fracture, / the distance between the inlet and outlet also called fracture length, and h the aperture of the fracture. 
The cubic law is a simplification of the Darcy's law in Eq. The comparison of the cubic and Darcy's law shows 
that the permeability of the fracture can be identified as ^ = h^/l2. This means that only the aperture is needed to 
calculate the permeability of a single fracture. 

5.7. Visco-capillary two-phase flow in fractures 

Visco-capillary two-phase flow occurs in a fracture network containing for example gas and water with a meniscus 
separating the two phases. It is described as a transport process whereby the water in the fracture or pore volume of 
a rock formation is displaced by gas under the influence of viscous and capillary forces |[35]| . Assuming again an 
idealized fracture network regarded as a network of tiny connected capillary tubes, for each tube Young-Laplace's 
equation is used to describe the governing physics of a capillary pressure /?cap forming across the interface between 
the two static fluids, gas and water, due to the surface tension 

2o-cos^ 

Pc^p = -J . (11) 
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cr denotes the surface or interfacial tension (for gas/water: cTg^ ^ 0.0073 Nm"^), 6 is the wetting angle, and h is the 
aperture value of the fracture. In other words, Young-Laplace represents the difference between gas and water pressure 
needed to displace the water from the initially fully saturated fracture network. Relating the capillary pressure and the 
water saturation in the fracture network gives the capillary pressure - water saturation relationship, commonly called 
water retention curve. 



5.2. Flow simulation in the channel model 

Several approaches have been used to solve the flow equations stated above. Considering the case of the DFN with 
a large number of fractures in 3D space, a simplification is inevitable. Mostly, the description in the 3D network is 
replaced by a capillary model (tubes) as stated before. One approach of a capillary model is called channel model and 
was proposed by Cacas et al. |20|. It holds the idea that the intersection of two fractures is schematized by a channel 
(tube) that joins their centers with an eff'ective hydraulic conductivity resulting from simple geometric arguments. It 
is based on the observation that in real fracture planes, flow does not occur over the entire surface of the fracture, but 
only in a limited number of channels within the fracture plane. It further assumes that each intersection of the fracture 
planes is connected to a fictitious node at the center of each circular fracture by a single equivalent straight channel 
starting from the center of the intersection. The exact geometry of the real channel is thereby omitted. In this work 
we use the centers of each intersection as fictitious nodes. The center node is left out, since the flow in the fracture 
will always follow the direct path from node to node. 

To solve the flow problem, a global system matrix for the flow fracture network, containing the elements (channels) 
and nodes from the channel model is constructed. A similar equation for the cubic law in Eq. ([TO]) can be expressed as 

^=^mP, (12) 

where 2 g denotes the volumetric flow rate vector, P the pressure field vector, y the dynamic viscosity and 
M G M'^x" the global system matrix. The permeability of the leakage path, seen as the permeability of an "equivalent" 
single fracture, can be identified as 

^path = , (13) 

where ^path is the equivalent or total aperture of the leakage path. To summarize, results using the channel model are 
a very simple approximation of flow through geological fracture networks. Nevertheless, the essential physics from 
the single fracture scale to the much larger fracture network scale is included in the described model in a consistent 
and also computationally feasible way. Note that the flow calculations are only performed on the backbones, that have 
to be found via the TMIP, as shown in the following. 



6. Isotropic vs. anisotropic fracture networks 

To address the efl'ect of anisotropy, trapping MIP (Ci = 1) on the specific nearly isotropic artificial 3D DFN 
(see Fig. [2] a)) and the highly anisotropic DFN (see Fig. [2]b)) is studied. The DFNs are similar with respect to the 
local fracture porosity and fracture area, only that the fracture orientation is changed by setting the Fisher dispersion 
parameter from a: = 1 for the isotropic case to a: = 50 for the highly anisotropic one. In all scenarios, the entry zone 
is defined as the xz-plane on the left side at j = 0, and the exit zone correspondingly as the xz-plane on the right side 
8ity = L. Periodic boundaries are applied on the other four sides. Accordingly, the invasion process starts at the entry 
zone and progresses in the direction of the outlet. Instead of stopping the invasion process after reaching the outlet, 
leading to the breakthrough in Fig. [5] the process continues by assuming an infinite gas production rate at the inlet. 
Finally, the process stops after full invasion in Fig.|6] 

The concentration in terms of the average number of intersections per fracture is calculated numerically as p = 
(Mnt) ~ 6.1 for the isotropic (a: = 1) and p ^ 1.93 for the highly anisotropic (a: = 50) DFN. Hence, the concentration 
of the anisotropic DFN is actually considerably smaller than the one of the isotropic DFN. The asymptotic values (for 
the infinite system size L ^ oo) found by Khamforoush et al. ||27]| for anisotropic DFNs with a: = 50 are p^oo^^ ^2.1 
and p[.ooy^ ~ 2.1 in x- and j-direction and p^co^^ ~ 2.4 in z-direction. Given that the MIP proceeds in j-direction 
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suggests a concentration value slightly below the asymptotic value. Finite size effects are expected and are discussed 
in Seed 

The considerably smaller amount of fracture intersections (concentration) is immediately apparent comparing the 
amount of invaded fractures after full invasion in Fig.|6]a) and d). A large part, most notably in the middle of the DFN, 
is not invaded, leaving a very high proportion amount of trapped water in the DFN (Fig. [6]c) and f)). The detailed 
description of the invasion process is investigated in Sec.ISj together with a real-world application. 
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Figure 5: DFN at breakthrough. 

a) and d) Invaded fractures marked from early (red) to late invasion (blue) after trapping MIP from left _y = to right y = L.b) and e) Flow fracture 
network (backbone of the percolation cluster) and channel model discretization of the non-pruned fractures, c) and f) Pressure distribution in the 
flow fracture network. The fractures are colored according to the mean pressure on the fracture plane. Red indicates high and blue low pressure, 
a), b), and c) correspond to an isotropic DFN (k = 1) while d), e), and f) are for an anisotropic DFN (k = 50). 
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Figure 6: DFN after full invasion with trapping MIR 

a) and d) Invaded fractures, b) and e) Pressure distribution in the flow fracture network, c) and f) Invaded fractures that contain trapped water. 
Colors indicate the amount of trapped water from fully trapped with 100% (dark blue) water content to non-trapped (fully invaded) with 0% (white) 
water content. Transparent fractures at the exit zone (representing the outlet) are not identified (no trapping possible). 



7. Finite size dependence 

To study size effects, trapping MIP is applied on several artificial 3D discrete fracture networks (DFN). The DFNs 
are generated in a unit cell of size L^. The polygons, representing the fractures, are oriented using a varying Fisher 
dispersion parameter k = {1, 10, 20, 30,40, 50}. Four diff'erent system sizes are considered with system lengths L = 
{5 <imax, 10 <imax, 20 <imax, 30 <imax}, with the diameter of the circumscribed circle of the largest fracture in the DFN (imax- 
To assure that the DFN models have a connectivity independent of L and between entry and exit, a fracture density of 
Pfr = 4 was chosen, leading to DFN models slightly above the percolation thresholds Pc of ordinary percolation (OP) 
for the chosen system sizes (L < 30 Jmax)- The amount of fractures A^fr for L = {5 <imax, 10 <imax, 20 <imax, 30 Jmax} 
are therefore Nf^ = {500,4000, 32000, 108000} given that Nfr = Pfr/L^ However, it turns out that the concentration 
in terms of the average number of intersections per fracture (p' = ^Mnt)) is in the range of the critical concentration 

for high K values. The mean average concentration ^p (/c)^ for diff'erent system sizes and varying k is calculated 
numerically and plotted in Fig. [7| Obviously, the resulting mean average concentrations for a: = 40 and a: = 50 are 
below the critical concentration p^ (illustrated by the the dashed line in Fig.j?]). Additionally, the water saturation after 
full invasion with trapping MIP of DFNs with = 10 t/max containing 4000 fractures is considered. The gray coloring 
in Fig. [7] shows a decreasing invaded or gas volume for increasing k value. The curve resembles the decreasing 
concentration curve for L = 10 Jmax satisfactorily. 
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Figure 7: Mean average concentration for DFNs with varying system sizes, averaged over 1000 (L\), 200 (L2), and 100 (L3 and L4) realizations. 



To study the influence of the Fisher dispersion parameter k on the dimensionless critical aperture threshold and 
on the dimensionless permeability of the flow path or backbone K'^^^^ as well as to determine the scaling exponents of 
the difl'erent system sizes, the finite- size dependence is analyzed. The fracture network permeability is expected to 
scale as |[36i 

K^ip^^L) - L-'l\ (14) 

with conductivity exponent t and correlation length exponent v. At the criticality, the exponents are believed to be 
universal, i.e., to depend only on the dimension of the system, but not on the underlying lattice properties |37|. For 
DFNs, Koudina et al. 136 1 report tlv = 2.22 ± 0.08, which is in good agreement with the value for 3D t/v ^ 2.26 | 38 1. 
Above the critical concentration, the exponent is no longer universal and we found for a: < 30 

4ath(^'^)~i"^^''^^ (15) 

with dimensionless path permeability K'^^^^ dependent on k and system length L. The dependency is given by the 
nonlinear scaling function sk(k). Similarly, a scaling law can be stated for the dimensionless critical aperture threshold 
with 

V>,L)~L-^^^>, (16) 

with the nonlinear scaling function Sh(K). The critical aperture threshold and the path permeability for the artificial 
DFN rescaled accordingly are in Figs.[8]a) and c), and the functions sk(k) and Sh(K) are in Figs.[8]b) and d). Trapping 
MIP without inclination adjustment shows a linear decrease for the path permeability and a nonlinear decrease for the 
critical aperture threshold. This is in agreement with the results for ordinary percolation (OP) found by Khamforoush 
etal |27|. 

Trapping MIP with inclination adjustment shows an extremely nonlinear behavior. The reason for this is the strong 
eff'ect of the inclination adjustment for fractures, that are nearly perpendicular to each other for a small k value. In 
other words, the inclination adjustment anticipates a higher critical aperture threshold and path permeability with a 
decreasing strength for an increasing k. This is the counterpart to reduction of the critical aperture threshold and 
path permeability with decreasing average concentration p (increasing a:) shown in Fig. [t] Both eff'ects combined 
constitute this concave function with a maximum between k = 20 and k = 30. 
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Figure 8: a) The critical aperture threshold /zjj^ ^ rescaled by U^^*^^ where ShiK) is in b). Trapping MIP with varying inclination adjustment constants 
Ci, one with no path length adjustment (np = 1), and ordinary trapping IP with np = nc = nj = 1 (small plot) are studied, c) Similarly, the path 
permeability ^p^^j^ at breakthrough on the same DFN rescaled by U^^'^^ with sk{j<) in d). 



8. Comparison to "real" fracture network 

The trapping MIP (d = 1) was applied to the DFN model of the specific Wellenberg realization with side length 
L = 100 m containing 12684 fractures. The invasion process started at the left hand site at j = m and progressed to 
the right hand side dXy - 100 m. The breakthrough and the full invasion state are shown in Fig. [9] The full invasion 
state of this specific Wellenberg realization is studied and compared to the artificially generated DFN models. 
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Figure 9: Wellenberg DFN at breakthrough and after full invasion. 

a) Invaded fractures at breakthrough marked from early (red) to late invasion (blue) after trapping MIP from left _y = to right _y = 100 m. b) Flow 
fracture network (backbone of the percolation cluster) and channel model discretization of the non-pruned fractures, c) Pressure distribution in the 
flow fracture network. The fractures are colored according to the mean pressure on the fracture plane. Red indicates high and blue low pressure, d) 
Invaded fractures after continuing the invading process till full invasion with trapping MIP. e) Pressure distribution in the flow fracture network, f) 
Invaded fractures that contain trapped water. Colors indicate the amount of trapped water from fully trapped with 100% (dark blue) water content 
to non-trapped (fully invaded) with 0% (white) water content. Transparent fractures at the exit zone (representing the outlet) are not identified (no 
trapping possible). 

The decline of the aperture threshold for each invasion step can be used to study the dynamics of the invading 
process. The decreasing aperture threshold illustrates thereby the rising entry pressure. Additionally, the sum of 
invasions into fractures at each invasion step called invasion rate describes the advance of the gas migration in the 
DFN. Note that a single fracture might have several inlets, represented by the connections to other fractures. Hence, 
more than one invasion into a single fracture is possible and the total sum of invasion rates typically exceeds the total 
amount of fractures for well connected DFNs. The aperture threshold value and the invasion rate for each step are 
shown in Fig.|9] Jmax is set to 1 m for the artificially generated DFNs, to be able to directly compare the results to the 
Wellenberg DFN. 

The breakthrough for the nearly isotropic DFN (a: = 1) is reached after 918 steps and at a critical aperture (invasion 
percolation) threshold value of K^^^ ^ 7.4110"^ m. The aperture threshold curve, plotted in logarithmic scale, has a 
constant slope of Z? ^ -1.27 after breakthrough till around 5000 steps. For the highly anisotropic DFN (k = 50), 
breakthrough is reached after 553 steps and at a critical aperture (invasion percolation) threshold value of h[^^ ^ 
4.7710"^ m. Hence the critical aperture threshold is more than one order of magnitude below the one of the isotropic 
case. After breakthrough only 25 more invasion steps are performed till the full invasion state is reached, where no 
more invasion is possible. This is caused by the lack of connectivity of the anisotropic DFN and emphasizes, that the 
system with infinite size is not percolating for < p^^. 

Apparently, a huge number of invasion steps have to be performed for the Wellenberg DFN till full invasion, due to 
its high connectivity. The invasion rate is defined as the sum of all invasions into fractures. The faults in the Wellenberg 
DFN model are tessellated into equally sized fractures. Hence, multiple invasions into fractures belonging to the same 
macro fracture (same fault) can occur at the same step. The natural behavior of taking the path of the least resistance 
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might lead to preferred spreading in macro fractures for a certain amount of steps. A closer inspection of Fig. [9] a) 
and b) with respect to the invasion step coloring of macro fractures supports this assumption. The resulting capillary 
pressure is calculated using Eq. ( 1J_) and plotted in relation to the water saturation. As assumed, only a small amount 
of trapped water remains in the Wellenberg DFN. The aperture threshold and water retention curves of the isotropic 
and the Wellenberg DFN are very similar, with slightly steeper constant slope (b ^ 1.71) between breakthrough and 
around 45000 steps for the Wellenberg DFN. Although the Wellenberg DFN is of anisotropic nature. Also the water 
retention curve (capillary pressure - water saturation curve) resembles the one of the nearly isotropic artificial DFN 
in both shape and values. The only diff'erence is that the remaining water saturation of the Wellenberg DFN at full 
invasion is about 10% smaller than in the artificial one. This confirms the above-mentioned assumption, that the 
higher connectivity of the DFN leads to a smaller amount of trapped water after full invasion. 
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Figure 10: On the left: The invasion process with the aperture threshold and the invasion rate for each step for trapping MIP on the three different 
DFNs. The breakthrough step is marked by a dashed line for the invasion steps and by a dotted line for the aperture threshold value. On the right: 
Capillary pressure - water saturation relationship with the total amount of trapped water marked by a dashed line. The amount of trapped water of 
pruned "dead end" fractures, that are not invadable (determined before the MIP) is indicated by a solid cross line (perpendicular to the saturation 
axis). 



9. Conclusions 

This work has investigated two-phase flow in anisotropic three-dimensional (3D) discrete fracture networks (DFN). 
The DFNs are either based on a general distribution function or have been conditioned on geological measurements 
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from boreholes, core and outcrop. A modified invasion percolation (MIP) was proposed to study percolation of the 
immiscible fluids through the DFN. The MIP was designed to describe the behavior at fracture intersections in more 
detail by including the contact angles and the length of the line of intersection between connected fractures. Further- 
more, it incorporates the hydraulic path length from one connection point to the other. 

A systematic study of the critical aperture threshold and the path permeability of the remaining flow fracture net- 
work after trapping MIP on artificial anisotropic 3D DFNs is shown in this work. The DFNs are constructed according 
to the Fisher distribution. The anisotropy is thereby controlled by the Fisher dispersion parameter k. With increasing 
/c, the distribution of the fracture normal vector becomes more dense and inclined towards the initial mean direction, 
chosen to be the z-axis. It is found for trapping MIP without inclination adjustment, that increased anisotropy in 
z-direction leads to a decreased critical aperture (invasion percolation) threshold. The dependency of the spanning 
path (infinite cluster) to the anisotropy is in good agreement with the study of Hunt et al. (391 and Khamforoush 
et al. ill] . Furthermore, the path permeability of the remaining flow fracture networks is decreased similarly with 
increased anisotropy. Including the inclination adjustment anticipates a higher critical aperture threshold and path per- 
meability with a decreasing strength for increasing k. This is the counterpart to the reduction of the critical aperture 
threshold and path permeability with increasing k. Both efl'ects combined constitute a concave function, that increases 
with K until a certain value depending on the average concentration in terms of the average number of intersections 
per fracture and decreases afterwards. Additional modifications of the invasion percolation procedure, like gravity 
adjustment, might have a similar influence on the critical aperture threshold and path permeability to k curves. These 
findings are in good agreement with our physical perception and can be used for extremely anisotropic DFNs. 

The "real" application, the so-called Wellenberg DFN model, was found to be directly comparable to the artificial 
3D DFN. Due to the high connectivity of the Wellenberg DFN model, it resembles mostly the nearly isotropic artificial 
DFN. The main diff'erence is found to be the remaining water saturation of the Wellenberg DFN at full invasion, which 
is smaller than in the artificial one. This is due to the higher connectivity of the DFN that avoids trapping during the 
MIP process. It was shown in Ref. | lOJ, that trapping in general 3D lattices is very rare. Trapping MIP on fully 
isotropic DFN with extremely high connectivity is expected to show the same characteristics. This was demonstrated 
for many applications of porous media 131] |371 |40l . Note that the isolated boundaries of the Wellenberg model in 
comparison to the periodic boundaries of the artificial DFNs are unproblematic due to the high connectivity of the 
Wellenberg. 

It is our belief, that invasion percolation with physical modifications is a useful numerical tool for the physical 
abstraction of the complicated situation of two-phase flow in fractured rock. It helps to reduce the problem to flow 
in a percolation backbone containing a very small portion of the total amount of fractures, that could potentially 
be invaded. Further studies should focus on the comparison with more detailed physics-based models, volumetric 
calculations with fluid dynamics of flow in typical situations like fracture intersections with measured geometries to 
advance the inclination adjustments, or two-phase flow inside real fracture cavities in specific rock as such. Since 
percolation features are strongly aff'ected by correlations in the media fTT-TTl, it would also be relevant to analyze the 
impact of spatially correlated distributions of apertures and fracture orientations on the flow properties in successive 
work. 
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